library(here)
## here() starts at /Users/julienvollering/Desktop/PortableOffice/conifer-plantation-spread
library(tidyverse)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.1     ✔ purrr   0.3.2
## ✔ tibble  2.1.1     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(glmmTMB)
library(DHARMa)
## Registered S3 method overwritten by 'DHARMa':
##   method        from   
##   refit.glmmTMB glmmTMB
datlist <- readRDS(here("output","data.rds"))
locality.preds <- c('age','bio01','bio19','seeds.ExP','seeds.WALD','relelev')
contrast <- "T4"
excluded <- "plantation forest"

P.sitchensis-lutz

Ps_dat <- Filter(function(x)!all(is.na(x)), datlist$`P.sitchensis-lutz`)
Ps_arith <- pivot_longer(Ps_dat, 10:ncol(Ps_dat), names_to = "nin", values_to = "are") %>%
  group_by(nin) %>% 
  summarize(daa = sum(are)/10, # units from are to decare
            n = sum(wildlings*are)) %>% 
  mutate(density = n/daa) %>%
  arrange(desc(density), desc(n), daa)
scaler <- filter(Ps_arith, nin == contrast) %>% 
  pull(density)
Ps_arith <- mutate(Ps_arith, fit = density/scaler)
nin daa n density fit
T17 2.0 42.8 21.1 5.4
T2 160.2 961.2 6.0 1.5
T34 1071.0 5664.4 5.3 1.3
T4 1362.9 5360.0 3.9 1.0
fen 758.3 2886.0 3.8 1.0
disturbed 409.5 1550.0 3.8 1.0
plantation forest 88.9 227.5 2.6 0.7
T32 676.5 1236.0 1.8 0.5
T16 14.8 25.3 1.7 0.4
swamp 160.5 251.0 1.6 0.4
bog 186.8 252.8 1.4 0.3
T6 34.9 39.4 1.1 0.3
T1 141.5 158.0 1.1 0.3
T27 68.1 75.3 1.1 0.3
T31 317.4 342.3 1.1 0.3
T21 6.3 6.4 1.0 0.3
T29 6.7 4.8 0.7 0.2
cultivated 1830.1 1015.8 0.6 0.1
T24 6.9 3.1 0.5 0.1
T12 11.7 4.6 0.4 0.1
T13 11.7 2.0 0.2 0.0
T18 2.0 0.0 0.0 0.0
T3 2.3 0.0 0.0 0.0
T33 5.9 0.0 0.0 0.0
spring 9.8 0.0 0.0 0.0
perfectlyseparated <- c('T18', 'T3', 'T33', 'spring')
removed <- c(perfectlyseparated, excluded)
Ps_dat <- Ps_dat %>% 
  filter_at(removed, all_vars(. < 0.5)) %>% 
  filter(!(is.na(seeds.ExP) | is.na(seeds.ExP))) %>% 
  select(-c('species'))
Ps_means <- Ps_dat[, locality.preds] %>% map_dbl(mean)
Ps_sds <- Ps_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
  Ps_dat <- modify_at(Ps_dat, .at = i, ~ (. - Ps_means[i])/Ps_sds[i])
}
summary(Ps_dat)
##    locality           wildlings             age              bio01        
##  Length:73309       Min.   :  0.0000   Min.   :-1.5919   Min.   :-2.3687  
##  Class :character   1st Qu.:  0.0000   1st Qu.:-0.9296   1st Qu.:-0.5921  
##  Mode  :character   Median :  0.0000   Median : 0.1006   Median :-0.3229  
##                     Mean   :  0.2788   Mean   : 0.0000   Mean   : 0.0000  
##                     3rd Qu.:  0.0000   3rd Qu.: 0.5422   3rd Qu.: 0.9082  
##                     Max.   :206.0000   Max.   : 2.2347   Max.   : 1.5957  
##      bio19           seeds.ExP         seeds.WALD          relelev       
##  Min.   :-1.8113   Min.   :-0.2245   Min.   :-0.49584   Min.   :-3.3307  
##  1st Qu.:-0.5378   1st Qu.:-0.2245   1st Qu.:-0.47446   1st Qu.:-0.4699  
##  Median :-0.2429   Median :-0.2245   Median :-0.40001   Median : 0.0833  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.9368   3rd Qu.:-0.2241   3rd Qu.:-0.05296   3rd Qu.: 0.4406  
##  Max.   : 2.5052   Max.   :13.1032   Max.   :14.58650   Max.   : 7.9503  
##  plantation forest     cultivated       disturbed            T32         
##  Min.   :0.0000000   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:0.0000000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000  
##  Median :0.0000000   Median :0.0000   Median :0.00000   Median :0.00000  
##  Mean   :0.0006327   Mean   :0.2495   Mean   :0.05579   Mean   :0.09212  
##  3rd Qu.:0.0000000   3rd Qu.:0.5250   3rd Qu.:0.00000   3rd Qu.:0.00000  
##  Max.   :0.4910000   Max.   :1.0000   Max.   :1.00000   Max.   :1.00000  
##       fen              T34              T13                 T4        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.000000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.000000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :0.000000   Median :0.0000  
##  Mean   :0.1033   Mean   :0.1458   Mean   :0.001598   Mean   :0.1858  
##  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.000000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.000000   Max.   :1.0000  
##        T2                T1              swamp             T18          
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.00e+00  
##  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.00e+00  
##  Median :0.00000   Median :0.00000   Median :0.0000   Median :0.00e+00  
##  Mean   :0.02185   Mean   :0.01929   Mean   :0.0219   Mean   :7.64e-05  
##  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.00e+00  
##  Max.   :1.00000   Max.   :1.00000   Max.   :1.0000   Max.   :4.67e-01  
##       T27                T17                spring        
##  Min.   :0.000000   Min.   :0.0000000   Min.   :0.000000  
##  1st Qu.:0.000000   1st Qu.:0.0000000   1st Qu.:0.000000  
##  Median :0.000000   Median :0.0000000   Median :0.000000  
##  Mean   :0.009282   Mean   :0.0002765   Mean   :0.001233  
##  3rd Qu.:0.000000   3rd Qu.:0.0000000   3rd Qu.:0.000000  
##  Max.   :0.999000   Max.   :1.0000000   Max.   :0.499000  
##       bog                T6                T31               T16          
##  Min.   :0.00000   Min.   :0.000000   Min.   :0.00000   Min.   :0.000000  
##  1st Qu.:0.00000   1st Qu.:0.000000   1st Qu.:0.00000   1st Qu.:0.000000  
##  Median :0.00000   Median :0.000000   Median :0.00000   Median :0.000000  
##  Mean   :0.02548   Mean   :0.004763   Mean   :0.04324   Mean   :0.002025  
##  3rd Qu.:0.00000   3rd Qu.:0.000000   3rd Qu.:0.00000   3rd Qu.:0.000000  
##  Max.   :1.00000   Max.   :1.000000   Max.   :1.00000   Max.   :1.000000  
##       T29                T12                T24          
##  Min.   :0.000000   Min.   :0.000000   Min.   :0.000000  
##  1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.000000  
##  Median :0.000000   Median :0.000000   Median :0.000000  
##  Mean   :0.000916   Mean   :0.001592   Mean   :0.000945  
##  3rd Qu.:0.000000   3rd Qu.:0.000000   3rd Qu.:0.000000  
##  Max.   :0.922000   Max.   :1.000000   Max.   :0.500000  
##        T3                T21                 T33          
##  Min.   :0.00e+00   Min.   :0.0000000   Min.   :0.000000  
##  1st Qu.:0.00e+00   1st Qu.:0.0000000   1st Qu.:0.000000  
##  Median :0.00e+00   Median :0.0000000   Median :0.000000  
##  Mean   :4.52e-05   Mean   :0.0008594   Mean   :0.000216  
##  3rd Qu.:0.00e+00   3rd Qu.:0.0000000   3rd Qu.:0.000000  
##  Max.   :4.27e-01   Max.   :1.0000000   Max.   :0.498000
paste(names(Ps_dat)[!(names(Ps_dat) %in% c(contrast, removed))], collapse = " + ")
## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + cultivated + disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 + T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21"

Compare WALD vs ExP

Ps_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + 
                 cultivated + disturbed + T32 + fen + T34 + T13 + T2 + T1 + 
                 swamp + T27 + T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + 
                 T21 + (1 | locality)")
Ps_m1 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ 0)
summary(Ps_m1)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +  
##     T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |  
##     locality)
## Data: Ps_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##  44876.3  45133.9 -22410.1  44820.3    73281 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 3.745    1.935   
## Number of obs: 73309, groups:  locality, 36
## 
## Overdispersion parameter for nbinom2 family (): 0.0788 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -2.3588408  0.3312308   -7.12 1.07e-12 ***
## seeds.WALD      0.9384915  0.0266383   35.23  < 2e-16 ***
## age             0.2583634  0.3360831    0.77 0.442043    
## bio01           0.5949349  0.4141349    1.44 0.150839    
## bio19          -0.5308042  0.4021168   -1.32 0.186827    
## relelev        -0.2617019  0.0360642   -7.26 3.97e-13 ***
## cultivated     -2.6215100  0.1013181  -25.87  < 2e-16 ***
## disturbed      -0.1931504  0.1223341   -1.58 0.114364    
## T32            -0.9567007  0.1188110   -8.05 8.13e-16 ***
## fen            -0.4023827  0.1093915   -3.68 0.000235 ***
## T34            -0.1604257  0.1097421   -1.46 0.143784    
## T13             0.3284157  0.8961519    0.37 0.714012    
## T2              0.0008741  0.1998395    0.00 0.996510    
## T1             -2.0711687  0.2804754   -7.38 1.53e-13 ***
## swamp          -2.3173703  0.2126976  -10.90  < 2e-16 ***
## T27             1.1609938  0.4601492    2.52 0.011633 *  
## T17             2.8532089  1.0404665    2.74 0.006102 ** 
## bog            -1.5504993  0.1971566   -7.86 3.71e-15 ***
## T6             -1.4496908  0.3818381   -3.80 0.000147 ***
## T31             1.0036357  0.2054121    4.89 1.03e-06 ***
## T16             0.3939334  0.4827710    0.82 0.414509    
## T29            -1.4845290  1.1672067   -1.27 0.203421    
## T12            -1.8898514  1.4458260   -1.31 0.191176    
## T24            -3.1043266  2.1623452   -1.44 0.151108    
## T21            -3.4083792  0.7636041   -4.46 8.06e-06 ***
## seeds.WALD:age  0.1599230  0.0259332    6.17 6.97e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ps_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + 
                 cultivated + disturbed + T32 + fen + T34 + T13 + T2 + T1 + 
                 swamp + T27 + T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + 
                 T21 + (1 | locality)")
Ps_m2 <- glmmTMB(Ps_f2, Ps_dat, family = "nbinom2", ziformula = ~ 0)
summary(Ps_m2)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +  
##     T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |  
##     locality)
## Data: Ps_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##  46062.6  46320.2 -23003.3  46006.6    73281 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 3.268    1.808   
## Number of obs: 73309, groups:  locality, 36
## 
## Overdispersion parameter for nbinom2 family (): 0.0639 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.00098    0.31012  -6.452 1.10e-10 ***
## seeds.ExP      0.52240    0.02614  19.985  < 2e-16 ***
## age            0.40449    0.31396   1.288 0.197629    
## bio01          0.50619    0.38742   1.307 0.191359    
## bio19         -0.59032    0.37635  -1.569 0.116761    
## relelev       -0.31592    0.03810  -8.292  < 2e-16 ***
## cultivated    -2.90729    0.10156 -28.625  < 2e-16 ***
## disturbed     -0.12715    0.12970  -0.980 0.326947    
## T32           -0.98601    0.11678  -8.444  < 2e-16 ***
## fen           -0.61490    0.11376  -5.405 6.47e-08 ***
## T34           -0.15135    0.11279  -1.342 0.179644    
## T13            0.14798    0.92685   0.160 0.873153    
## T2             0.29751    0.21278   1.398 0.162057    
## T1            -2.77607    0.28435  -9.763  < 2e-16 ***
## swamp         -2.62667    0.22223 -11.819  < 2e-16 ***
## T27            0.74253    0.50100   1.482 0.138316    
## T17            3.22023    1.19005   2.706 0.006811 ** 
## bog           -1.34171    0.19957  -6.723 1.78e-11 ***
## T6            -1.30451    0.40234  -3.242 0.001186 ** 
## T31            1.05021    0.20777   5.055 4.31e-07 ***
## T16            0.43948    0.51149   0.859 0.390216    
## T29           -0.81335    1.20018  -0.678 0.497965    
## T12           -2.57557    1.44967  -1.777 0.075624 .  
## T24           -0.67897    1.88438  -0.360 0.718614    
## T21           -2.71096    0.76065  -3.564 0.000365 ***
## seeds.ExP:age -0.03212    0.02600  -1.235 0.216686    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##       dAIC   df
## Ps_m1    0.0 28
## Ps_m2 1186.3 28

Continue with WALD.

Compare zero inflation effects

Ps_m3 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ 1)
summary(Ps_m3)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +  
##     T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |  
##     locality)
## Zero inflation:             ~1
## Data: Ps_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##  44878.3  45145.1 -22410.1  44820.3    73280 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 3.745    1.935   
## Number of obs: 73309, groups:  locality, 36
## 
## Overdispersion parameter for nbinom2 family (): 0.0788 
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -2.3588826  0.3312219   -7.12 1.07e-12 ***
## seeds.WALD      0.9385074  0.0266387   35.23  < 2e-16 ***
## age             0.2582882  0.3360732    0.77 0.442162    
## bio01           0.5948216  0.4141232    1.44 0.150906    
## bio19          -0.5308007  0.4021053   -1.32 0.186817    
## relelev        -0.2616964  0.0360642   -7.26 3.98e-13 ***
## cultivated     -2.6215078  0.1013185  -25.87  < 2e-16 ***
## disturbed      -0.1932046  0.1223341   -1.58 0.114263    
## T32            -0.9566731  0.1188114   -8.05 8.14e-16 ***
## fen            -0.4023916  0.1093918   -3.68 0.000235 ***
## T34            -0.1604635  0.1097424   -1.46 0.143691    
## T13             0.3284436  0.8961246    0.37 0.713980    
## T2              0.0008862  0.1998404    0.00 0.996462    
## T1             -2.0711986  0.2804756   -7.38 1.53e-13 ***
## swamp          -2.3173310  0.2126978  -10.89  < 2e-16 ***
## T27             1.1607030  0.4601308    2.52 0.011651 *  
## T17             2.8532586  1.0404906    2.74 0.006102 ** 
## bog            -1.5504769  0.1971573   -7.86 3.72e-15 ***
## T6             -1.4496469  0.3818395   -3.80 0.000147 ***
## T31             1.0035677  0.2054108    4.89 1.03e-06 ***
## T16             0.3940391  0.4827768    0.82 0.414390    
## T29            -1.4845403  1.1672111   -1.27 0.203419    
## T12            -1.8896773  1.4457986   -1.31 0.191208    
## T24            -3.1045989  2.1623385   -1.44 0.151071    
## T21            -3.4085395  0.7636025   -4.46 8.05e-06 ***
## seeds.WALD:age  0.1599326  0.0259335    6.17 6.96e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -17.33     550.26  -0.031    0.975
Ps_m4 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ age)
summary(Ps_m4)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +  
##     T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |  
##     locality)
## Zero inflation:             ~age
## Data: Ps_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##  44355.8  44631.9 -22147.9  44295.8    73279 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 3.744    1.935   
## Number of obs: 73309, groups:  locality, 36
## 
## Overdispersion parameter for nbinom2 family (): 0.16 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.51775    0.33423   -4.54 5.60e-06 ***
## seeds.WALD      0.90437    0.02792   32.40  < 2e-16 ***
## age            -0.25650    0.33716   -0.76 0.446792    
## bio01           0.61522    0.41432    1.48 0.137576    
## bio19          -0.53728    0.40203   -1.34 0.181411    
## relelev        -0.26289    0.03411   -7.71 1.28e-14 ***
## cultivated     -2.58291    0.10399  -24.84  < 2e-16 ***
## disturbed      -0.43691    0.12949   -3.37 0.000741 ***
## T32            -1.09632    0.11718   -9.36  < 2e-16 ***
## fen            -0.39970    0.11135   -3.59 0.000331 ***
## T34            -0.14437    0.10765   -1.34 0.179893    
## T13             0.32869    0.89216    0.37 0.712565    
## T2             -0.87796    0.24094   -3.64 0.000269 ***
## T1             -2.34639    0.31287   -7.50 6.40e-14 ***
## swamp          -2.38736    0.25218   -9.47  < 2e-16 ***
## T27             1.42264    0.48746    2.92 0.003517 ** 
## T17             2.50901    0.77584    3.23 0.001221 ** 
## bog            -1.70283    0.22673   -7.51 5.90e-14 ***
## T6             -2.03677    0.39124   -5.21 1.93e-07 ***
## T31             1.03591    0.20327    5.10 3.47e-07 ***
## T16             0.40353    0.45464    0.89 0.374768    
## T29            -2.16010    1.48667   -1.45 0.146228    
## T12            -2.35583    1.83304   -1.29 0.198723    
## T24            -3.34399    2.79879   -1.19 0.232165    
## T21            -3.28725    0.72425   -4.54 5.66e-06 ***
## seeds.WALD:age  0.11690    0.02579    4.53 5.84e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.14693    0.09179  -1.601    0.109    
## age         -1.11316    0.05366 -20.747   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ps_m5 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(Ps_m5)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +  
##     T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |  
##     locality)
## Zero inflation:             ~age + (1 | locality)
## Data: Ps_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##  42266.6  42551.9 -21102.3  42204.6    73278 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 2.034    1.426   
## Number of obs: 73309, groups:  locality, 36
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 5.133    2.266   
## Number of obs: 73309, groups:  locality, 36
## 
## Overdispersion parameter for nbinom2 family (): 0.274 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -0.35045    0.26151   -1.34 0.180208    
## seeds.WALD      0.83919    0.02646   31.72  < 2e-16 ***
## age             0.04685    0.26979    0.17 0.862133    
## bio01          -0.15971    0.40960   -0.39 0.696587    
## bio19          -0.86385    0.38061   -2.27 0.023230 *  
## relelev        -0.30122    0.03738   -8.06 7.72e-16 ***
## cultivated     -2.63220    0.10448  -25.19  < 2e-16 ***
## disturbed      -0.76252    0.15240   -5.00 5.63e-07 ***
## T32            -0.80398    0.11738   -6.85 7.40e-12 ***
## fen            -0.15901    0.11135   -1.43 0.153271    
## T34             0.03014    0.10599    0.28 0.776104    
## T13             0.47073    0.91836    0.51 0.608249    
## T2              1.04641    0.26711    3.92 8.95e-05 ***
## T1             -2.32908    0.24152   -9.64  < 2e-16 ***
## swamp          -2.02971    0.29586   -6.86 6.86e-12 ***
## T27             1.56716    0.42690    3.67 0.000242 ***
## T17             2.30316    0.63314    3.64 0.000275 ***
## bog            -1.35952    0.27339   -4.97 6.60e-07 ***
## T6              0.72195    0.74465    0.97 0.332290    
## T31             0.73905    0.23772    3.11 0.001877 ** 
## T16             0.51663    0.50709    1.02 0.308294    
## T29            -2.27767    1.85349   -1.23 0.219127    
## T12            -2.79530    1.87336   -1.49 0.135664    
## T24            -1.86526    2.58922   -0.72 0.471283    
## T21            -5.10913    1.31469   -3.89 0.000102 ***
## seeds.WALD:age  0.18578    0.02372    7.83 4.78e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   1.2857     0.3924   3.277  0.00105 **
## age          -0.8856     0.3880  -2.282  0.02247 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##       dAIC   df
## Ps_m5    0.0 31
## Ps_m4 2089.2 30
## Ps_m1 2609.6 28
## Ps_m3 2611.6 29

m5 is better than m1, m3, and m4: much better AIC, clear zero-inflation effects

simulationOutput <- simulateResiduals(fittedModel = Ps_m5)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.075717, p-value = 0.104
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0202, p-value = 0.248
## alternative hypothesis: two.sided

Diagnostics look OK.

Interpret relative establishment probability

selmod <- Ps_m5
terms <- attributes(selmod$frame)$terms 
newdat <- attributes(terms)$factors %>% 
  as_tibble() %>% 
  select(-`seeds.WALD:age`) %>% 
  mutate(seeds.WALD = (Ps_means["seeds.WALD"] - Ps_means["seeds.WALD"]) / Ps_sds["seeds.WALD"],
         age = (Ps_means["age"] - Ps_means["age"]) / Ps_sds["age"],
         bio01 = (Ps_means["bio01"] - Ps_means["bio01"]) / Ps_sds["bio01"],
         bio19 = (Ps_means["bio19"] - Ps_means["bio19"]) / Ps_sds["bio19"],
         relelev = (Ps_means["relelev"] - Ps_means["relelev"]) / Ps_sds["relelev"],
         locality = NA_character_) %>% 
  distinct()
newdat <- newdat %>% 
  mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>% 
  replace_na(list(nin = "T4"))
Ps_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
Ps_link <- mutate(Ps_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>% 
  select(-se.fit) %>% 
  mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(Ps_link, nin == "T4") %>% 
  pull(fit)
link.scaled <- mutate_at(Ps_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))
gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) + 
  geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) + 
  geom_point()  + 
  geom_point(data = filter(Ps_arith, nin %in% link.scaled$nin), 
             col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()

newdat.locality <- newdat %>% 
  select(-locality) %>% 
  expand_grid(locality = c(NA_character_, unique(Ps_dat$locality)))
resp.locality <- newdat.locality %>% 
  bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>% 
  pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))

gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) + 
  geom_point(data = filter(resp.locality.scaled, !is.na(locality)), 
             mapping = aes(color = locality), cex = 0.5, show.legend = FALSE)  + 
  geom_point(data = filter(resp.locality.scaled, is.na(locality)))  + 
  geom_point(data = filter(Ps_arith, nin %in% resp.locality.scaled$nin), 
             color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()

Larix

L_dat <- Filter(function(x)!all(is.na(x)), datlist$Larix)
L_arith <- pivot_longer(L_dat, 10:ncol(L_dat), names_to = "nin", values_to = "are") %>%
  group_by(nin) %>% 
  summarize(daa = sum(are)/10, # units from are to decare
            n = sum(wildlings*are)) %>% 
  mutate(density = n/daa) %>%
  arrange(desc(density), desc(n), daa)
scaler <- filter(L_arith, nin == contrast) %>% 
  pull(density)
L_arith <- mutate(L_arith, fit = density/scaler)
nin daa n density fit
T17 2.0 211.8 104.5 160.6
plantation forest 169.1 159.8 0.9 1.5
T34 471.2 308.7 0.7 1.0
T4 332.4 216.3 0.7 1.0
T32 227.5 115.0 0.5 0.8
T1 40.3 6.9 0.2 0.3
disturbed 185.8 17.5 0.1 0.1
cultivated 1243.6 91.1 0.1 0.1
fen 132.8 6.1 0.0 0.1
T18 0.6 0.0 0.0 0.0
T30 1.0 0.0 0.0 0.0
T27 1.5 0.0 0.0 0.0
T2 1.6 0.0 0.0 0.0
swamp 1.7 0.0 0.0 0.0
T13 3.0 0.0 0.0 0.0
perfectlyseparated <- c('T18', 'T30', 'T2', 'T27', 'swamp', 'T13')
removed <- c(perfectlyseparated, excluded)
L_dat <- L_dat %>% 
  filter_at(removed, all_vars(. < 0.5)) %>% 
  filter(!(is.na(seeds.ExP) | is.na(seeds.ExP))) %>% 
  select(-c('species'))
L_means <- L_dat[, locality.preds] %>% map_dbl(mean)
L_sds <- L_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
  L_dat <- modify_at(L_dat, .at = i, ~ (. - L_means[i])/L_sds[i])
}
summary(L_dat)
##    locality           wildlings             age         
##  Length:26514       Min.   : 0.00000   Min.   :-1.9249  
##  Class :character   1st Qu.: 0.00000   1st Qu.:-0.2941  
##  Mode  :character   Median : 0.00000   Median : 0.3353  
##                     Mean   : 0.03719   Mean   : 0.0000  
##                     3rd Qu.: 0.00000   3rd Qu.: 0.5642  
##                     Max.   :27.00000   Max.   : 1.5655  
##      bio01              bio19           seeds.ExP         seeds.WALD     
##  Min.   :-2.53170   Min.   :-1.5815   Min.   :-0.1189   Min.   :-0.3842  
##  1st Qu.: 0.05949   1st Qu.:-0.4381   1st Qu.:-0.1189   1st Qu.:-0.3666  
##  Median : 0.29676   Median :-0.2804   Median :-0.1189   Median :-0.3109  
##  Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.67912   3rd Qu.: 0.4687   3rd Qu.:-0.1189   3rd Qu.:-0.1421  
##  Max.   : 0.90390   Max.   : 2.2034   Max.   :21.3893   Max.   :19.1019  
##     relelev         plantation forest   cultivated       disturbed      
##  Min.   :-3.16949   Min.   :0.00000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:-0.54556   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Median :-0.03968   Median :0.00000   Median :0.1405   Median :0.00000  
##  Mean   : 0.00000   Mean   :0.00384   Mean   :0.4678   Mean   :0.06954  
##  3rd Qu.: 0.46001   3rd Qu.:0.00000   3rd Qu.:1.0000   3rd Qu.:0.00000  
##  Max.   : 5.73614   Max.   :0.49900   Max.   :1.0000   Max.   :1.00000  
##       T32               fen               T34              T13           
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000000  
##  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000000  
##  Median :0.00000   Median :0.00000   Median :0.0000   Median :0.0000000  
##  Mean   :0.08533   Mean   :0.04989   Mean   :0.1766   Mean   :0.0002235  
##  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000000  
##  Max.   :1.00000   Max.   :1.00000   Max.   :1.0000   Max.   :0.4650000  
##        T4               T2                  T1              swamp         
##  Min.   :0.0000   Min.   :0.0000000   Min.   :0.00000   Min.   :0.00e+00  
##  1st Qu.:0.0000   1st Qu.:0.0000000   1st Qu.:0.00000   1st Qu.:0.00e+00  
##  Median :0.0000   Median :0.0000000   Median :0.00000   Median :0.00e+00  
##  Mean   :0.1247   Mean   :0.0002907   Mean   :0.01493   Mean   :6.74e-05  
##  3rd Qu.:0.0000   3rd Qu.:0.0000000   3rd Qu.:0.00000   3rd Qu.:0.00e+00  
##  Max.   :1.0000   Max.   :0.4980000   Max.   :1.00000   Max.   :4.83e-01  
##       T30         T18         T27                T17           
##  Min.   :0   Min.   :0   Min.   :0.000000   Min.   :0.0000000  
##  1st Qu.:0   1st Qu.:0   1st Qu.:0.000000   1st Qu.:0.0000000  
##  Median :0   Median :0   Median :0.000000   Median :0.0000000  
##  Mean   :0   Mean   :0   Mean   :0.000159   Mean   :0.0007646  
##  3rd Qu.:0   3rd Qu.:0   3rd Qu.:0.000000   3rd Qu.:0.0000000  
##  Max.   :0   Max.   :0   Max.   :0.419000   Max.   :1.0000000
paste(names(L_dat)[!(names(L_dat) %in% c(contrast, removed))], collapse = " + ")
## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + cultivated + disturbed + T32 + fen + T34 + T1 + T17"

Compare WALD vs ExP

L_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + 
                cultivated + disturbed + T32 + fen + T34 + T1 + T17 + 
                (1 | locality)")
L_m1 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ 0)
summary(L_m1)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Data: L_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3878.1   4009.1  -1923.1   3846.1    26498 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.2738   0.5233  
## Number of obs: 26514, groups:  locality, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.0595 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -3.84395    0.29882 -12.864  < 2e-16 ***
## seeds.WALD      1.03433    0.06119  16.903  < 2e-16 ***
## age            -0.44714    0.31379  -1.425 0.154170    
## bio01          -1.03506    0.39701  -2.607 0.009130 ** 
## bio19          -1.93007    0.47899  -4.029 5.59e-05 ***
## relelev         0.02552    0.08020   0.318 0.750277    
## cultivated     -2.43001    0.32461  -7.486 7.10e-14 ***
## disturbed      -2.41780    0.46959  -5.149 2.62e-07 ***
## T32            -1.22154    0.32713  -3.734 0.000188 ***
## fen            -3.72529    0.83550  -4.459 8.24e-06 ***
## T34             0.03735    0.35197   0.106 0.915484    
## T1             -1.94085    0.97643  -1.988 0.046844 *  
## T17             6.40605    1.05478   6.073 1.25e-09 ***
## seeds.WALD:age  0.24491    0.06875   3.562 0.000367 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + 
                cultivated + disturbed + T32 + fen + T34 + T1 + T17 + 
                (1 | locality)")
L_m2 <- glmmTMB(L_f2, L_dat, family = "nbinom2", ziformula = ~ 0)
summary(L_m2)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Data: L_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   4249.0   4379.9  -2108.5   4217.0    26498 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.7066   0.8406  
## Number of obs: 26514, groups:  locality, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.0331 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.42176    0.35463  -6.829 8.55e-12 ***
## seeds.ExP      0.64949    0.08092   8.026 1.01e-15 ***
## age            0.41651    0.43886   0.949 0.342584    
## bio01          0.21247    0.56587   0.375 0.707303    
## bio19         -0.56967    0.67811  -0.840 0.400861    
## relelev       -0.03232    0.09079  -0.356 0.721811    
## cultivated    -3.63714    0.34925 -10.414  < 2e-16 ***
## disturbed     -3.27412    0.46437  -7.051 1.78e-12 ***
## T32           -1.72551    0.34102  -5.060 4.20e-07 ***
## fen           -4.50569    0.65557  -6.873 6.29e-12 ***
## T34           -1.36713    0.37758  -3.621 0.000294 ***
## T1            -1.92174    0.82430  -2.331 0.019734 *  
## T17            5.04949    1.33503   3.782 0.000155 ***
## seeds.ExP:age -0.07588    0.08216  -0.924 0.355696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##      dAIC  df
## L_m1   0.0 16
## L_m2 370.8 16

Continue with WALD.

Compare zero inflation effects

L_m3 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ 1)
summary(L_m3)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Zero inflation:             ~1
## Data: L_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3880.1   4019.3  -1923.1   3846.1    26497 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.2738   0.5233  
## Number of obs: 26514, groups:  locality, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.0595 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -3.84392    0.29882 -12.864  < 2e-16 ***
## seeds.WALD      1.03432    0.06119  16.903  < 2e-16 ***
## age            -0.44713    0.31379  -1.425 0.154181    
## bio01          -1.03507    0.39702  -2.607 0.009130 ** 
## bio19          -1.93010    0.47900  -4.029 5.59e-05 ***
## relelev         0.02552    0.08020   0.318 0.750276    
## cultivated     -2.43007    0.32461  -7.486 7.10e-14 ***
## disturbed      -2.41785    0.46959  -5.149 2.62e-07 ***
## T32            -1.22158    0.32713  -3.734 0.000188 ***
## fen            -3.72532    0.83550  -4.459 8.24e-06 ***
## T34             0.03731    0.35197   0.106 0.915583    
## T1             -1.94086    0.97642  -1.988 0.046843 *  
## T17             6.40607    1.05479   6.073 1.25e-09 ***
## seeds.WALD:age  0.24490    0.06875   3.562 0.000367 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -16.7     2484.9  -0.007    0.995
L_m4 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ age)
summary(L_m4)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Zero inflation:             ~age
## Data: L_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3869.9   4017.3  -1917.0   3833.9    26496 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.2046   0.4523  
## Number of obs: 26514, groups:  locality, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.0644 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -3.52811    0.30007 -11.758  < 2e-16 ***
## seeds.WALD      1.02067    0.06336  16.109  < 2e-16 ***
## age            -0.98463    0.32757  -3.006  0.00265 ** 
## bio01          -1.19452    0.36823  -3.244  0.00118 ** 
## bio19          -2.06985    0.43459  -4.763 1.91e-06 ***
## relelev         0.06530    0.09487   0.688  0.49125    
## cultivated     -2.42125    0.31614  -7.659 1.88e-14 ***
## disturbed      -2.43422    0.46514  -5.233 1.66e-07 ***
## T32            -1.26467    0.32281  -3.918 8.94e-05 ***
## fen            -3.79127    0.82779  -4.580 4.65e-06 ***
## T34            -0.02954    0.34640  -0.085  0.93204    
## T1             -1.86342    0.94799  -1.966  0.04934 *  
## T17             6.39796    1.01722   6.290 3.18e-10 ***
## seeds.WALD:age  0.25390    0.09144   2.777  0.00549 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -3.400      2.975  -1.143    0.253
## age           -2.458      1.537  -1.599    0.110
L_m5 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(L_m5)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Zero inflation:             ~age + (1 | locality)
## Data: L_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3843.2   3998.7  -1902.6   3805.2    26495 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.2603   0.5102  
## Number of obs: 26514, groups:  locality, 12
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 1.411    1.188   
## Number of obs: 26514, groups:  locality, 12
## 
## Overdispersion parameter for nbinom2 family (): 0.138 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -2.77051    0.36671  -7.555 4.19e-14 ***
## seeds.WALD      1.04741    0.06627  15.804  < 2e-16 ***
## age            -1.23721    0.42699  -2.898  0.00376 ** 
## bio01          -0.96633    0.43191  -2.237  0.02526 *  
## bio19          -2.07193    0.52678  -3.933 8.38e-05 ***
## relelev         0.04995    0.09175   0.544  0.58612    
## cultivated     -2.45366    0.32924  -7.452 9.17e-14 ***
## disturbed      -2.39472    0.48467  -4.941 7.78e-07 ***
## T32            -1.37959    0.33335  -4.139 3.49e-05 ***
## fen            -3.53032    0.80225  -4.401 1.08e-05 ***
## T34             0.14793    0.36357   0.407  0.68409    
## T1             -2.02332    0.87769  -2.305  0.02115 *  
## T17             6.05552    0.71886   8.424  < 2e-16 ***
## seeds.WALD:age  0.26766    0.10735   2.493  0.01265 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.4100     0.5420   0.756   0.4494  
## age          -1.0206     0.4953  -2.061   0.0393 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##      dAIC df
## L_m5  0.0 19
## L_m4 26.7 18
## L_m1 34.9 16
## L_m3 36.9 17

m5 is better than m1, m3, and m4: much better AIC, clear zero-inflation effects

simulationOutput <- simulateResiduals(fittedModel = L_m5)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.25346, p-value = 0.144
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99942, p-value = 0.816
## alternative hypothesis: two.sided

Diagnostics look OK.

Interpret relative establishment probability

selmod <- L_m5
terms <- attributes(selmod$frame)$terms 
newdat <- attributes(terms)$factors %>% 
  as_tibble() %>% 
  select(-`seeds.WALD:age`) %>% 
  mutate(seeds.WALD = (Ps_means["seeds.WALD"] - L_means["seeds.WALD"]) / L_sds["seeds.WALD"],
         age = (Ps_means["age"] - L_means["age"]) / L_sds["age"],
         bio01 = (Ps_means["bio01"] - L_means["bio01"]) / L_sds["bio01"],
         bio19 = (Ps_means["bio19"] - L_means["bio19"]) / L_sds["bio19"],
         relelev = (Ps_means["relelev"] - L_means["relelev"]) / L_sds["relelev"],
         locality = NA_character_) %>% 
  distinct()
newdat <- newdat %>% 
  mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>% 
  replace_na(list(nin = "T4"))
L_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
L_link <- mutate(L_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>% 
  select(-se.fit) %>% 
  mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(L_link, nin == "T4") %>% 
  pull(fit)
link.scaled <- mutate_at(L_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))
gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) + 
  geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) + 
  geom_point()  + 
  geom_point(data = filter(L_arith, nin %in% link.scaled$nin), 
             col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()

newdat.locality <- newdat %>% 
  select(-locality) %>% 
  expand_grid(locality = c(NA_character_, unique(L_dat$locality)))
resp.locality <- newdat.locality %>% 
  bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>% 
  pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))

gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) + 
  geom_point(data = filter(resp.locality.scaled, !is.na(locality)), 
             mapping = aes(color = locality), cex = 0.5, show.legend = FALSE)  + 
  geom_point(data = filter(resp.locality.scaled, is.na(locality)))  + 
  geom_point(data = filter(L_arith, nin %in% resp.locality.scaled$nin), 
             color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()

P.abies

Pa_dat <- Filter(function(x)!all(is.na(x)), datlist$P.abies)
Pa_arith <- pivot_longer(Pa_dat, 10:ncol(Pa_dat), names_to = "nin", values_to = "are") %>%
  group_by(nin) %>% 
  summarize(daa = sum(are)/10, # units from are to decare
            n = sum(wildlings*are)) %>% 
  mutate(density = n/daa) %>%
  arrange(desc(density), desc(n), daa)
scaler <- filter(Pa_arith, nin == contrast) %>% 
  pull(density)
Pa_arith <- mutate(Pa_arith, fit = density/scaler)
nin daa n density fit
T1 13.8 34.0 2.5 2.0
T4 689.3 843.4 1.2 1.0
swamp 33.2 28.6 0.9 0.7
plantation forest 14.5 10.5 0.7 0.6
T32 306.8 49.9 0.2 0.1
fen 96.8 14.6 0.2 0.1
spring 6.1 0.8 0.1 0.1
disturbed 76.0 7.9 0.1 0.1
bog 11.8 1.0 0.1 0.1
cultivated 631.3 4.9 0.0 0.0
T2 17.7 0.1 0.0 0.0
T6 2.9 0.0 0.0 0.0
T31 17.7 0.0 0.0 0.0
perfectlyseparated <- c('T6', 'T31', 'T2')
removed <- c(perfectlyseparated, excluded)
Pa_dat <- Pa_dat %>% 
  filter_at(removed, all_vars(. < 0.5)) %>% 
  filter(!(is.na(seeds.ExP) | is.na(seeds.ExP))) %>% 
  select(-c('species'))
Pa_means <- Pa_dat[, locality.preds] %>% map_dbl(mean)
Pa_sds <- Pa_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
  Pa_dat <- modify_at(Pa_dat, .at = i, ~ (. - Pa_means[i])/Pa_sds[i])
}
summary(Pa_dat)
##    locality           wildlings            age              bio01        
##  Length:18786       Min.   : 0.0000   Min.   :-1.4707   Min.   :-2.4383  
##  Class :character   1st Qu.: 0.0000   1st Qu.:-0.7900   1st Qu.:-0.4425  
##  Mode  :character   Median : 0.0000   Median :-0.2999   Median : 0.3878  
##                     Mean   : 0.0528   Mean   : 0.0000   Mean   : 0.0000  
##                     3rd Qu.: 0.0000   3rd Qu.: 0.4352   3rd Qu.: 0.6571  
##                     Max.   :13.0000   Max.   : 1.8238   Max.   : 0.8414  
##      bio19           seeds.ExP         seeds.WALD       
##  Min.   :-1.7911   Min.   :-0.2017   Min.   :-0.579549  
##  1st Qu.:-0.3616   1st Qu.:-0.2017   1st Qu.:-0.497256  
##  Median : 0.1501   Median :-0.2017   Median :-0.371481  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.000000  
##  3rd Qu.: 1.0679   3rd Qu.:-0.2014   3rd Qu.: 0.002906  
##  Max.   : 1.2141   Max.   :14.5311   Max.   :10.142177  
##     relelev         plantation forest     cultivated       disturbed      
##  Min.   :-4.65908   Min.   :0.0000000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:-0.35233   1st Qu.:0.0000000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Median : 0.04963   Median :0.0000000   Median :0.0000   Median :0.00000  
##  Mean   : 0.00000   Mean   :0.0003784   Mean   :0.3357   Mean   :0.04047  
##  3rd Qu.: 0.51376   3rd Qu.:0.0000000   3rd Qu.:1.0000   3rd Qu.:0.00000  
##  Max.   : 3.39527   Max.   :0.4850000   Max.   :1.0000   Max.   :1.00000  
##       T32              fen                T4               T2          
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.000000  
##  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.000000  
##  Median :0.0000   Median :0.00000   Median :0.0000   Median :0.000000  
##  Mean   :0.1632   Mean   :0.05135   Mean   :0.3661   Mean   :0.002631  
##  3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:1.0000   3rd Qu.:0.000000  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :0.499000  
##        T1               swamp            spring              bog          
##  Min.   :0.000000   Min.   :0.0000   Min.   :0.000000   Min.   :0.000000  
##  1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.:0.000000   1st Qu.:0.000000  
##  Median :0.000000   Median :0.0000   Median :0.000000   Median :0.000000  
##  Mean   :0.003123   Mean   :0.0177   Mean   :0.003221   Mean   :0.005946  
##  3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.:0.000000   3rd Qu.:0.000000  
##  Max.   :1.000000   Max.   :1.0000   Max.   :1.000000   Max.   :1.000000  
##        T6                 T31           
##  Min.   :0.0000000   Min.   :0.0000000  
##  1st Qu.:0.0000000   1st Qu.:0.0000000  
##  Median :0.0000000   Median :0.0000000  
##  Mean   :0.0005545   Mean   :0.0004149  
##  3rd Qu.:0.0000000   3rd Qu.:0.0000000  
##  Max.   :0.4980000   Max.   :0.4480000
paste(names(Pa_dat)[!(names(Pa_dat) %in% c(contrast, removed))], collapse = " + ")
## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + cultivated + disturbed + T32 + fen + T1 + swamp + spring + bog"

Compare WALD vs ExP

Pa_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + 
                 cultivated + disturbed + T32 + fen + T1 + swamp + spring + bog + 
                 (1 | locality)")
Pa_m1 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pa_m1)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   4160.2   4293.5  -2063.1   4126.2    18769 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 1.827    1.352   
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for nbinom2 family (): 0.347 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -3.28264    0.46979  -6.987 2.80e-12 ***
## seeds.WALD      0.44241    0.04472   9.894  < 2e-16 ***
## age            -0.42684    0.75175  -0.568    0.570    
## bio01           0.23004    0.60744   0.379    0.705    
## bio19           0.70574    0.86238   0.818    0.413    
## relelev        -0.32952    0.06404  -5.145 2.67e-07 ***
## cultivated     -6.66521    0.58588 -11.376  < 2e-16 ***
## disturbed      -2.88120    0.46898  -6.144 8.07e-10 ***
## T32            -1.39290    0.19489  -7.147 8.86e-13 ***
## fen            -1.73769    0.42175  -4.120 3.78e-05 ***
## T1             -0.15064    0.42250  -0.357    0.721    
## swamp          -1.29798    0.32585  -3.983 6.79e-05 ***
## spring         -2.22515    1.54449  -1.441    0.150    
## bog            -1.09333    1.37612  -0.795    0.427    
## seeds.WALD:age  0.05965    0.05143   1.160    0.246    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pa_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + 
                 cultivated + disturbed + T32 + fen + T1 + swamp + spring + bog +
                 (1 | locality)")
Pa_m2 <- glmmTMB(Pa_f2, Pa_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pa_m2)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   4275.2   4408.5  -2120.6   4241.2    18769 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 1.861    1.364   
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for nbinom2 family (): 0.347 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.08756    0.47254  -6.534 6.41e-11 ***
## seeds.ExP     -0.01568    0.06770  -0.232    0.817    
## age           -0.41037    0.75759  -0.542    0.588    
## bio01          0.16883    0.61227   0.276    0.783    
## bio19          0.79135    0.86941   0.910    0.363    
## relelev       -0.36921    0.06050  -6.103 1.04e-09 ***
## cultivated    -6.90142    0.60046 -11.494  < 2e-16 ***
## disturbed     -2.98567    0.44660  -6.685 2.30e-11 ***
## T32           -1.49012    0.19098  -7.802 6.07e-15 ***
## fen           -1.80618    0.41803  -4.321 1.56e-05 ***
## T1            -0.26371    0.41785  -0.631    0.528    
## swamp         -1.42100    0.32366  -4.390 1.13e-05 ***
## spring        -1.77133    1.47632  -1.200    0.230    
## bog           -1.26941    1.42697  -0.890    0.374    
## seeds.ExP:age -0.09854    0.09317  -1.058    0.290    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##       dAIC df
## Pa_m1   0  17
## Pa_m2 115  17

Continue with WALD.

Compare zero inflation effects

Pa_m3 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ 1)
summary(Pa_m3)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Zero inflation:             ~1
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   4162.2   4303.4  -2063.1   4126.2    18768 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 1.827    1.352   
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for nbinom2 family (): 0.347 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -3.28263    0.46979  -6.987 2.80e-12 ***
## seeds.WALD      0.44241    0.04472   9.894  < 2e-16 ***
## age            -0.42684    0.75175  -0.568    0.570    
## bio01           0.23004    0.60744   0.379    0.705    
## bio19           0.70574    0.86238   0.818    0.413    
## relelev        -0.32952    0.06404  -5.145 2.67e-07 ***
## cultivated     -6.66522    0.58588 -11.376  < 2e-16 ***
## disturbed      -2.88121    0.46898  -6.144 8.07e-10 ***
## T32            -1.39290    0.19489  -7.147 8.86e-13 ***
## fen            -1.73770    0.42175  -4.120 3.78e-05 ***
## T1             -0.15064    0.42250  -0.357    0.721    
## swamp          -1.29798    0.32585  -3.983 6.79e-05 ***
## spring         -2.22504    1.54444  -1.441    0.150    
## bog            -1.09325    1.37609  -0.794    0.427    
## seeds.WALD:age  0.05965    0.05143   1.160    0.246    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -18.88    2199.10  -0.009    0.993
Pa_m4 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ age)
summary(Pa_m4)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Zero inflation:             ~age
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   4082.0   4230.9  -2022.0   4044.0    18767 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 1.645    1.283   
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for nbinom2 family (): 0.563 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.70065    0.49367  -3.445 0.000571 ***
## seeds.WALD      0.47071    0.08818   5.338 9.40e-08 ***
## age             1.05118    0.75423   1.394 0.163404    
## bio01           0.14233    0.58405   0.244 0.807471    
## bio19           0.53330    0.82029   0.650 0.515608    
## relelev        -0.28260    0.06231  -4.535 5.75e-06 ***
## cultivated     -6.66195    0.58627 -11.363  < 2e-16 ***
## disturbed      -2.97487    0.47862  -6.216 5.12e-10 ***
## T32            -1.39752    0.19406  -7.201 5.96e-13 ***
## fen            -1.87243    0.40895  -4.579 4.68e-06 ***
## T1             -0.12570    0.38750  -0.324 0.745634    
## swamp          -1.22709    0.30478  -4.026 5.67e-05 ***
## spring         -2.36531    1.56487  -1.512 0.130660    
## bog            -0.99254    1.36192  -0.729 0.466135    
## seeds.WALD:age  0.05025    0.10939   0.459 0.645939    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.6100     0.2981   2.046   0.0407 *  
## age           2.5112     0.2385  10.528   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pa_m5 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(Pa_m5)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Zero inflation:             ~age + (1 | locality)
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3787.2   3944.0  -1873.6   3747.2    18766 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.2872   0.5359  
## Number of obs: 18786, groups:  locality, 9
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 5.737    2.395   
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for nbinom2 family (): 0.974 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.16094    0.32224  -3.603 0.000315 ***
## seeds.WALD      0.66246    0.09936   6.668 2.60e-11 ***
## age            -0.37772    0.47387  -0.797 0.425401    
## bio01          -0.40545    0.33732  -1.202 0.229381    
## bio19           0.32013    0.44298   0.723 0.469877    
## relelev        -0.19049    0.07535  -2.528 0.011474 *  
## cultivated     -6.72680    0.58479 -11.503  < 2e-16 ***
## disturbed      -3.40811    0.52840  -6.450 1.12e-10 ***
## T32            -1.61926    0.19840  -8.162 3.31e-16 ***
## fen            -1.32823    0.39136  -3.394 0.000689 ***
## T1             -0.10743    0.33060  -0.325 0.745208    
## swamp          -1.10088    0.27419  -4.015 5.94e-05 ***
## spring         -3.90949    1.74967  -2.234 0.025455 *  
## bog            -1.70301    1.37601  -1.238 0.215850    
## seeds.WALD:age  0.01866    0.12540   0.149 0.881723    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   1.6569     0.8732   1.898   0.0578 .
## age           1.0916     0.8401   1.299   0.1938  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##       dAIC  df
## Pa_m5   0.0 20
## Pa_m4 294.8 19
## Pa_m1 373.0 17
## Pa_m3 375.0 18

m5 is better than m1, m3, and m4: much better AIC, zero-inflation effects

Drop insignificant interaction term?

Pa_f6 <- update.formula(Pa_f1, ~ . - seeds.WALD:age)
Pa_m6 <- glmmTMB(Pa_f6, Pa_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(Pa_m6)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD + age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Zero inflation:             ~age + (1 | locality)
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3785.2   3934.2  -1873.6   3747.2    18767 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.2927   0.5411  
## Number of obs: 18786, groups:  locality, 9
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 5.711    2.39    
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for nbinom2 family (): 0.975 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.15211    0.31936  -3.608 0.000309 ***
## seeds.WALD   0.65391    0.08092   8.081 6.43e-16 ***
## age         -0.35620    0.45463  -0.783 0.433335    
## bio01       -0.40959    0.33850  -1.210 0.226275    
## bio19        0.32207    0.44547   0.723 0.469680    
## relelev     -0.19114    0.07519  -2.542 0.011021 *  
## cultivated  -6.73032    0.58456 -11.514  < 2e-16 ***
## disturbed   -3.42030    0.52325  -6.537 6.29e-11 ***
## T32         -1.61740    0.19790  -8.173 3.01e-16 ***
## fen         -1.32466    0.39052  -3.392 0.000694 ***
## T1          -0.10484    0.33013  -0.318 0.750801    
## swamp       -1.09958    0.27405  -4.012 6.01e-05 ***
## spring      -3.93216    1.74370  -2.255 0.024129 *  
## bog         -1.71428    1.37360  -1.248 0.212023    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   1.6642     0.8705   1.912   0.0559 .
## age           1.1061     0.8336   1.327   0.1845  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##       dAIC df
## Pa_m6  0   19
## Pa_m5  2   20
simulationOutput <- simulateResiduals(fittedModel = Pa_m6)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.050344, p-value = 0.016
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0206, p-value = 0.304
## alternative hypothesis: two.sided

Diagnostics look bad — particularly dispersion test.

Pa_m7 <- glmmTMB(Pa_f6, Pa_dat, family = "genpois", ziformula = ~ age + (1 | locality))
summary(Pa_m7)
##  Family: genpois  ( log )
## Formula:          
## wildlings ~ seeds.WALD + age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Zero inflation:             ~age + (1 | locality)
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3778.7   3927.7  -1870.4   3740.7    18767 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.09272  0.3045  
## Number of obs: 18786, groups:  locality, 9
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 7.579    2.753   
## Number of obs: 18786, groups:  locality, 9
## 
## Overdispersion parameter for genpois family (): 2.06 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.86769    0.24637  -7.581 3.44e-14 ***
## seeds.WALD   0.40904    0.04779   8.559  < 2e-16 ***
## age         -1.20711    0.37856  -3.189  0.00143 ** 
## bio01        0.35056    0.23679   1.480  0.13875    
## bio19        0.28991    0.36360   0.797  0.42526    
## relelev     -0.20621    0.07015  -2.940  0.00328 ** 
## cultivated  -6.30372    0.58887 -10.705  < 2e-16 ***
## disturbed   -4.04454    0.68253  -5.926 3.11e-09 ***
## T32         -1.47338    0.18570  -7.934 2.12e-15 ***
## fen         -1.30504    0.40327  -3.236  0.00121 ** 
## T1          -0.41010    0.35165  -1.166  0.24352    
## swamp       -0.98527    0.25745  -3.827  0.00013 ***
## spring      -1.95071    1.42452  -1.369  0.17088    
## bog         -1.12067    1.24805  -0.898  0.36922    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.01715    1.18071  -0.014    0.988
## age         -0.86046    1.11967  -0.768    0.442
simulationOutput <- simulateResiduals(fittedModel = Pa_m7)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.41878, p-value = 0.136
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0141, p-value = 0.464
## alternative hypothesis: two.sided

Better diagnostics with generalized poisson — including dispersion test.

Pa_m8 <- glmmTMB(Pa_f6, Pa_dat, family = poisson, ziformula = ~ age + (1 | locality))
summary(Pa_m8)
##  Family: poisson  ( log )
## Formula:          
## wildlings ~ seeds.WALD + age + bio01 + bio19 + relelev + cultivated +  
##     disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |  
##     locality)
## Zero inflation:             ~age + (1 | locality)
## Data: Pa_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3908.5   4049.6  -1936.2   3872.5    18768 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. 
##  locality (Intercept) 7.368e-10 2.714e-05
## Number of obs: 18786, groups:  locality, 9
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 3.031    1.741   
## Number of obs: 18786, groups:  locality, 9
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.008332   0.119825  -0.070 0.944562    
## seeds.WALD   0.218690   0.030478   7.175 7.21e-13 ***
## age         -0.255051   0.179221  -1.423 0.154706    
## bio01       -0.581891   0.126824  -4.588 4.47e-06 ***
## bio19        0.480586   0.160385   2.996 0.002731 ** 
## relelev      0.126163   0.056431   2.236 0.025372 *  
## cultivated  -6.693873   0.582473 -11.492  < 2e-16 ***
## disturbed   -2.850806   0.435732  -6.543 6.05e-11 ***
## T32         -1.728755   0.168677 -10.249  < 2e-16 ***
## fen         -1.068647   0.369595  -2.891 0.003835 ** 
## T1          -0.149395   0.288148  -0.518 0.604134    
## swamp       -0.920005   0.241456  -3.810 0.000139 ***
## spring      -2.619841   1.395613  -1.877 0.060491 .  
## bog         -1.659463   1.247526  -1.330 0.183451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.0831     0.6029   5.114 3.15e-07 ***
## age           0.8860     0.5953   1.488    0.137    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulationOutput <- simulateResiduals(fittedModel = Pa_m8)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.65524, p-value = 0.232
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.014, p-value = 0.616
## alternative hypothesis: two.sided

Better diagnostics with poisson — including dispersion test.

##       dAIC  df
## Pa_m7   0.0 19
## Pa_m6   6.5 19
## Pa_m8 129.8 18

Continue with generalized poisson model: best AIC and OK diagnostics.

Interpret relative establishment probability

selmod <- Pa_m7
terms <- attributes(selmod$frame)$terms 
newdat <- attributes(terms)$factors %>% 
  as_tibble() %>% 
  mutate(seeds.WALD = (Ps_means["seeds.WALD"] - Pa_means["seeds.WALD"]) / Pa_sds["seeds.WALD"],
         age = (Ps_means["age"] - Pa_means["age"]) / Pa_sds["age"],
         bio01 = (Ps_means["bio01"] - Pa_means["bio01"]) / Pa_sds["bio01"],
         bio19 = (Ps_means["bio19"] - Pa_means["bio19"]) / Pa_sds["bio19"],
         relelev = (Ps_means["relelev"] - Pa_means["relelev"]) / Pa_sds["relelev"],
         locality = NA_character_) %>% 
  distinct()
newdat <- newdat %>% 
  mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>% 
  replace_na(list(nin = "T4"))
Pa_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
Pa_link <- mutate(Pa_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>% 
  select(-se.fit) %>% 
  mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(Pa_link, nin == "T4") %>% 
  pull(fit)
link.scaled <- mutate_at(Pa_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))
gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) + 
  geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) + 
  geom_point()  + 
  geom_point(data = filter(Pa_arith, nin %in% link.scaled$nin), 
             col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()

newdat.locality <- newdat %>% 
  select(-locality) %>% 
  expand_grid(locality = c(NA_character_, unique(Pa_dat$locality)))
resp.locality <- newdat.locality %>% 
  bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>% 
  pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))

gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) + 
  geom_point(data = filter(resp.locality.scaled, !is.na(locality)), 
             mapping = aes(color = locality), cex = 0.5, show.legend = FALSE)  + 
  geom_point(data = filter(resp.locality.scaled, is.na(locality)))  + 
  geom_point(data = filter(Pa_arith, nin %in% resp.locality.scaled$nin), 
             color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()

P.contorta

Pc_dat <- Filter(function(x)!all(is.na(x)), datlist$P.contorta)
Pc_arith <- pivot_longer(Pc_dat, 10:ncol(Pc_dat), names_to = "nin", values_to = "are") %>%
  group_by(nin) %>% 
  summarize(daa = sum(are)/10, # units from are to decare
            n = sum(wildlings*are)) %>% 
  mutate(density = n/daa) %>%
  arrange(desc(density), desc(n), daa)
scaler <- filter(Pc_arith, nin == contrast) %>% 
  pull(density)
Pc_arith <- mutate(Pc_arith, fit = density/scaler)
nin daa n density fit
disturbed 216.2 3603.7 16.7 25.2
T34 117.7 233.6 2.0 3.0
T4 578.5 383.2 0.7 1.0
T30 2.5 1.6 0.6 1.0
T1 27.5 11.0 0.4 0.6
swamp 4.5 1.5 0.3 0.5
fen 219.1 43.4 0.2 0.3
plantation forest 119.1 12.1 0.1 0.2
bog 90.5 3.9 0.0 0.1
cultivated 7.6 0.0 0.0 0.0
T27 8.0 0.0 0.0 0.0
T32 23.3 0.0 0.0 0.0
perfectlyseparated <- c('cultivated', 'T27', 'T32')
removed <- c(perfectlyseparated, excluded)
Pc_dat <- Pc_dat %>% 
  filter_at(removed, all_vars(. < 0.5)) %>% 
  filter(!(is.na(seeds.ExP) | is.na(seeds.ExP) | is.na(relelev))) %>% 
  select(-c('species'))
Pc_means <- Pc_dat[, locality.preds] %>% map_dbl(mean)
Pc_sds <- Pc_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
  Pc_dat <- modify_at(Pc_dat, .at = i, ~ (. - Pc_means[i])/Pc_sds[i])
}
summary(Pc_dat)
##    locality           wildlings             age           
##  Length:12754       Min.   :  0.0000   Min.   :-1.413778  
##  Class :character   1st Qu.:  0.0000   1st Qu.:-0.539537  
##  Mode  :character   Median :  0.0000   Median : 0.006863  
##                     Mean   :  0.3473   Mean   : 0.000000  
##                     3rd Qu.:  0.0000   3rd Qu.: 0.334704  
##                     Max.   :131.0000   Max.   : 1.755345  
##      bio01             bio19           seeds.ExP         seeds.WALD     
##  Min.   :-1.0009   Min.   :-1.0620   Min.   :-0.1918   Min.   :-0.2757  
##  1st Qu.:-0.9670   1st Qu.:-1.0331   1st Qu.:-0.1918   1st Qu.:-0.2751  
##  Median :-0.8958   Median :-0.6807   Median :-0.1918   Median :-0.2701  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 1.1100   3rd Qu.: 0.8907   3rd Qu.:-0.1918   3rd Qu.:-0.2185  
##  Max.   : 1.1246   Max.   : 1.3528   Max.   :14.3222   Max.   :20.4881  
##     relelev        plantation forest    cultivated         disturbed     
##  Min.   :-2.2553   Min.   :0.000000   Min.   :0.000000   Min.   :0.0000  
##  1st Qu.:-0.4599   1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.0000  
##  Median :-0.1993   Median :0.000000   Median :0.000000   Median :0.0000  
##  Mean   : 0.0000   Mean   :0.004167   Mean   :0.000744   Mean   :0.1665  
##  3rd Qu.: 0.2092   3rd Qu.:0.000000   3rd Qu.:0.000000   3rd Qu.:0.0000  
##  Max.   : 5.1241   Max.   :0.499000   Max.   :0.496000   Max.   :1.0000  
##       T32                fen              T34                T4       
##  Min.   :0.000000   Min.   :0.0000   Min.   :0.00000   Min.   :0.000  
##  1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.000  
##  Median :0.000000   Median :0.0000   Median :0.00000   Median :0.014  
##  Mean   :0.001691   Mean   :0.1715   Mean   :0.09193   Mean   :0.451  
##  3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:1.000  
##  Max.   :0.495000   Max.   :1.0000   Max.   :1.00000   Max.   :1.000  
##        T1              swamp               T30          
##  Min.   :0.00000   Min.   :0.000000   Min.   :0.000000  
##  1st Qu.:0.00000   1st Qu.:0.000000   1st Qu.:0.000000  
##  Median :0.00000   Median :0.000000   Median :0.000000  
##  Mean   :0.02158   Mean   :0.003533   Mean   :0.001921  
##  3rd Qu.:0.00000   3rd Qu.:0.000000   3rd Qu.:0.000000  
##  Max.   :1.00000   Max.   :1.000000   Max.   :1.000000  
##       T27                 bog        
##  Min.   :0.0000000   Min.   :0.0000  
##  1st Qu.:0.0000000   1st Qu.:0.0000  
##  Median :0.0000000   Median :0.0000  
##  Mean   :0.0005656   Mean   :0.0707  
##  3rd Qu.:0.0000000   3rd Qu.:0.0000  
##  Max.   :0.4940000   Max.   :1.0000
paste(names(Pc_dat)[!(names(Pc_dat) %in% c(contrast, removed))], collapse = " + ")
## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + disturbed + fen + T34 + T1 + swamp + T30 + bog"

Compare WALD vs ExP

Pc_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + 
                 disturbed + fen + T34 + T1 + swamp + T30 + bog + (1 | locality)")
Pc_m1 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pc_m1)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + disturbed +  
##     fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Data: Pc_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3424.0   3543.3  -1696.0   3392.0    12738 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.06456  0.2541  
## Number of obs: 12754, groups:  locality, 6
## 
## Overdispersion parameter for nbinom2 family (): 0.0235 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -6.1381     0.4437 -13.834  < 2e-16 ***
## seeds.WALD       1.1642     0.1759   6.619 3.63e-11 ***
## age              0.1507     0.2331   0.647   0.5178    
## bio01           -8.5144     1.2909  -6.596 4.23e-11 ***
## bio19            6.3433     1.3481   4.705 2.53e-06 ***
## relelev         -0.1538     0.1986  -0.774   0.4388    
## disturbed        3.5411     0.3085  11.479  < 2e-16 ***
## fen              2.5660     1.0236   2.507   0.0122 *  
## T34              5.3311     0.9792   5.444 5.20e-08 ***
## T1               6.1689     1.0413   5.924 3.14e-09 ***
## swamp           -0.7505     2.5676  -0.292   0.7701    
## T30              1.4172     1.7734   0.799   0.4242    
## bog             -3.4078     0.7741  -4.402 1.07e-05 ***
## seeds.WALD:age  -0.5273     0.1211  -4.354 1.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pc_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + 
                 disturbed + fen + T34 + T1 + swamp + T30 + bog + (1 | locality)")
Pc_m2 <- glmmTMB(Pc_f2, Pc_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pc_m2)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + disturbed +  
##     fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Data: Pc_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3458.6   3577.9  -1713.3   3426.6    12738 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.04917  0.2217  
## Number of obs: 12754, groups:  locality, 6
## 
## Overdispersion parameter for nbinom2 family (): 0.022 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -6.026105   0.424305 -14.202  < 2e-16 ***
## seeds.ExP      0.801158   0.162047   4.944 7.65e-07 ***
## age            0.007757   0.216311   0.036 0.971393    
## bio01         -9.267952   1.195267  -7.754 8.91e-15 ***
## bio19          7.309543   1.259363   5.804 6.47e-09 ***
## relelev       -0.234139   0.197567  -1.185 0.235972    
## disturbed      3.650548   0.310504  11.757  < 2e-16 ***
## fen            2.376049   0.991193   2.397 0.016523 *  
## T34            5.458067   0.941805   5.795 6.82e-09 ***
## T1             6.447843   0.998741   6.456 1.08e-10 ***
## swamp          0.008496   2.252186   0.004 0.996990    
## T30            1.771512   1.802595   0.983 0.325727    
## bog           -3.361874   0.775087  -4.337 1.44e-05 ***
## seeds.ExP:age -0.357946   0.104856  -3.414 0.000641 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##       dAIC df
## Pc_m1  0.0 16
## Pc_m2 34.6 16

Continue with WALD.

Compare zero inflation effects

Pc_m3 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ 1)
summary(Pc_m3)
##  Family: nbinom2  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + disturbed +  
##     fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Zero inflation:             ~1
## Data: Pc_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3426.0   3552.7  -1696.0   3392.0    12737 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.06457  0.2541  
## Number of obs: 12754, groups:  locality, 6
## 
## Overdispersion parameter for nbinom2 family (): 0.0235 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -6.1381     0.4437 -13.834  < 2e-16 ***
## seeds.WALD       1.1642     0.1759   6.619 3.63e-11 ***
## age              0.1507     0.2331   0.647   0.5179    
## bio01           -8.5145     1.2909  -6.596 4.23e-11 ***
## bio19            6.3434     1.3481   4.705 2.53e-06 ***
## relelev         -0.1538     0.1986  -0.774   0.4388    
## disturbed        3.5411     0.3085  11.479  < 2e-16 ***
## fen              2.5660     1.0236   2.507   0.0122 *  
## T34              5.3311     0.9792   5.444 5.20e-08 ***
## T1               6.1689     1.0413   5.924 3.14e-09 ***
## swamp           -0.7506     2.5676  -0.292   0.7700    
## T30              1.4172     1.7734   0.799   0.4242    
## bog             -3.4078     0.7741  -4.402 1.07e-05 ***
## seeds.WALD:age  -0.5273     0.1211  -4.354 1.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -16.23    2437.67  -0.007    0.995
Pc_m4 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ age)
# Warning messages:
# 1: In fitTMB(TMBStruc) :
# Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
Pc_m5 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
# Warning messages:
# 1: In fitTMB(TMBStruc) :
# Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
# 2: In fitTMB(TMBStruc) :
# Model convergence problem; false convergence (8). See vignette('troubleshooting')

Overparameterized models? bio01 and bio19 show largest coefficients (5-6) while insignificant for other species. Try dropping them.

Pc_f6 <- update.formula(Pc_f1, ~ . - bio01 - bio19)
Pc_m6 <- glmmTMB(Pc_f6, Pc_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
# Warning messages:
# 1: In fitTMB(TMBStruc) :
# Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
Pc_m7 <- glmmTMB(Pc_f1, Pc_dat, family = "genpois", ziformula = ~ age + (1 | locality))
summary(Pc_m7)
##  Family: genpois  ( log )
## Formula:          
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + disturbed +  
##     fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Zero inflation:             ~age + (1 | locality)
## Data: Pc_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   3355.0   3496.6  -1658.5   3317.0    12735 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  locality (Intercept) 0.1903   0.4363  
## Number of obs: 12754, groups:  locality, 6
## 
## Zero-inflation model:
##  Groups   Name        Variance  Std.Dev. 
##  locality (Intercept) 6.884e-09 8.297e-05
## Number of obs: 12754, groups:  locality, 6
## 
## Overdispersion parameter for genpois family (): 87.5 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -4.12923    0.37565 -10.992  < 2e-16 ***
## seeds.WALD      0.30719    0.02676  11.480  < 2e-16 ***
## age             0.36613    0.29156   1.256   0.2092    
## bio01          -8.92063    1.57850  -5.651 1.59e-08 ***
## bio19           7.63388    1.53405   4.976 6.48e-07 ***
## relelev        -0.17693    0.11906  -1.486   0.1373    
## disturbed       3.38649    0.23107  14.655  < 2e-16 ***
## fen             0.71455    0.35603   2.007   0.0447 *  
## T34             3.52341    0.37936   9.288  < 2e-16 ***
## T1              5.49293    0.77380   7.099 1.26e-12 ***
## swamp           1.41629    1.62537   0.871   0.3836    
## T30             2.52351    1.03456   2.439   0.0147 *  
## bog            -1.66871    0.65506  -2.547   0.0109 *  
## seeds.WALD:age -0.03926    0.02091  -1.878   0.0604 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -4.161      8.316  -0.500    0.617
## age            2.955      4.737   0.624    0.533

Huge overdispersion parameter. Also huge coefficients for bio01 and bio19.

##       dAIC df
## Pc_m7  0   19
## Pc_m3 71   17
simulationOutput <- simulateResiduals(fittedModel = Pc_m7)
## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.83723, p-value = 0.648
## alternative hypothesis: two.sided
testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99957, p-value = 0.864
## alternative hypothesis: two.sided

Diagnostics look OK.

Continue with generalized poisson model: best AIC and OK diagnostics.

Interpret relative establishment probability

selmod <- Pc_m7
terms <- attributes(selmod$frame)$terms 
newdat <- attributes(terms)$factors %>% 
  as_tibble() %>% 
  select(-`seeds.WALD:age`) %>% 
  mutate(seeds.WALD = (Ps_means["seeds.WALD"] - Pc_means["seeds.WALD"]) / Pc_sds["seeds.WALD"],
         age = (Ps_means["age"] - Pc_means["age"]) / Pc_sds["age"],
         bio01 = (Ps_means["bio01"] - Pc_means["bio01"]) / Pc_sds["bio01"],
         bio19 = (Ps_means["bio19"] - Pc_means["bio19"]) / Pc_sds["bio19"],
         relelev = (Ps_means["relelev"] - Pc_means["relelev"]) / Pc_sds["relelev"],
         locality = NA_character_) %>% 
  distinct()
newdat <- newdat %>% 
  mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>% 
  replace_na(list(nin = "T4"))
Pc_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
Pc_link <- mutate(Pc_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>% 
  select(-se.fit) %>% 
  mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(Pc_link, nin == "T4") %>% 
  pull(fit)
link.scaled <- mutate_at(Pc_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))
gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) + 
  geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) + 
  geom_point()  + 
  geom_point(data = filter(Pc_arith, nin %in% link.scaled$nin), 
             col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()

newdat.locality <- newdat %>% 
  select(-locality) %>% 
  expand_grid(locality = c(NA_character_, unique(Pc_dat$locality)))
resp.locality <- newdat.locality %>% 
  bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>% 
  pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))

gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) + 
  geom_point(data = filter(resp.locality.scaled, !is.na(locality)), 
             mapping = aes(color = locality), cex = 0.5, show.legend = FALSE)  + 
  geom_point(data = filter(resp.locality.scaled, is.na(locality)))  + 
  geom_point(data = filter(L_arith, nin %in% resp.locality.scaled$nin), 
             color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()
## Warning: Transformation introduced infinite values in continuous y-axis

All species

link.list <- list(Ps_link, L_link, Pa_link, Pc_link)
names(link.list) <- c('Picea sitchensis/lutzii', 'Larix', 'Picea abies', 'Pinus contorta')
link <- bind_rows(link.list, .id = 'species')
scaler <- filter(link, species == 'Picea sitchensis/lutzii', nin == "T4") %>% 
  pull(fit)
link.scaled <- mutate_at(link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))
arith.list <- list(filter(Ps_arith, nin %in% Ps_link$nin), 
                   filter(L_arith, nin %in% L_link$nin),
                   filter(Pa_arith, nin %in% Pa_link$nin),
                   filter(Pc_arith, nin %in% Pc_link$nin))
names(arith.list) <- c('Picea sitchensis/lutzii', 'Larix', 'Picea abies', 'Pinus contorta')
arith <- bind_rows(arith.list, .id = 'species')
scaler <- filter(arith, species == 'Picea sitchensis/lutzii', nin == "T4") %>% 
  pull(density)
arith.scaled <- mutate(arith, fit = density/scaler)
link.scaled[,c('species','nin','fit','ci.l','ci.u')]
## # A tibble: 45 x 5
##    species                 nin           fit   ci.l   ci.u
##    <chr>                   <chr>       <dbl>  <dbl>  <dbl>
##  1 Picea sitchensis/lutzii T4         1      0.599   1.67 
##  2 Picea sitchensis/lutzii cultivated 0.0719 0.0432  0.120
##  3 Picea sitchensis/lutzii disturbed  0.466  0.268   0.811
##  4 Picea sitchensis/lutzii T32        0.448  0.264   0.760
##  5 Picea sitchensis/lutzii fen        0.853  0.510   1.43 
##  6 Picea sitchensis/lutzii T34        1.03   0.618   1.72 
##  7 Picea sitchensis/lutzii T13        1.60   0.250  10.3  
##  8 Picea sitchensis/lutzii T2         2.85   1.42    5.70 
##  9 Picea sitchensis/lutzii T1         0.0974 0.0502  0.189
## 10 Picea sitchensis/lutzii swamp      0.131  0.0625  0.276
## # … with 35 more rows
gg <- ggplot(link.scaled, aes(fct_reorder(nin, fit, first), fit, color = species)) + 
  geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u), 
                 position = position_dodge(width = 0.5)) + 
  geom_point(position = position_dodge(width = 0.5))  + 
  geom_point(data = arith.scaled, pch = 1, position = position_dodge(width = 0.5))
gg + scale_y_log10() + coord_flip() + 
  labs(x = 'NiN type', y = 'relative establishment probability', title = 'Relative establishment across species and NiN types')

Filled points are model estimates with lines showing 95% confidence intervals. Open points are calculated directly from the number of wildlings per unit area.

sessionInfo

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] DHARMa_0.2.4    glmmTMB_0.2.3   forcats_0.4.0   stringr_1.4.0  
##  [5] dplyr_0.8.3     purrr_0.3.2     readr_1.3.1     tidyr_1.0.0    
##  [9] tibble_2.1.1    ggplot2_3.1.1   tidyverse_1.2.1 here_0.1       
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.0          jsonlite_1.6        splines_3.6.0      
##  [4] foreach_1.4.4       modelr_0.1.4        gap_1.2.1          
##  [7] assertthat_0.2.1    highr_0.8           stats4_3.6.0       
## [10] cellranger_1.1.0    yaml_2.2.0          numDeriv_2016.8-1.1
## [13] pillar_1.4.0        backports_1.1.4     lattice_0.20-38    
## [16] glue_1.3.1          bbmle_1.0.22        digest_0.6.19      
## [19] rvest_0.3.4         minqa_1.2.4         colorspace_1.4-1   
## [22] htmltools_0.4.0     Matrix_1.2-17       plyr_1.8.4         
## [25] spaMM_3.0.0         pkgconfig_2.0.2     broom_0.5.2        
## [28] haven_2.1.0         mvtnorm_1.0-10      scales_1.0.0       
## [31] lme4_1.1-21         proxy_0.4-23        generics_0.0.2     
## [34] ellipsis_0.3.0      withr_2.1.2         pbapply_1.4-2      
## [37] TMB_1.7.15          lazyeval_0.2.2      cli_1.1.0          
## [40] magrittr_1.5        crayon_1.3.4        readxl_1.3.1       
## [43] evaluate_0.13       fansi_0.4.0         nlme_3.1-139       
## [46] MASS_7.3-51.4       xml2_1.2.0          tools_3.6.0        
## [49] hms_0.4.2           lifecycle_0.1.0     munsell_0.5.0      
## [52] compiler_3.6.0      rlang_0.4.2         grid_3.6.0         
## [55] nloptr_1.2.1        iterators_1.0.10    rstudioapi_0.10    
## [58] rmarkdown_1.12      boot_1.3-22         gtable_0.3.0       
## [61] codetools_0.2-16    R6_2.4.0            lubridate_1.7.4    
## [64] knitr_1.23          bdsmatrix_1.3-4     zeallot_0.1.0      
## [67] utf8_1.1.4          rprojroot_1.3-2     stringi_1.4.3      
## [70] parallel_3.6.0      Rcpp_1.0.1          vctrs_0.2.1        
## [73] tidyselect_0.2.5    xfun_0.7